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Resonant Scattering and Ly-alpha Radiation Emergent from 

Neutral Hydrogen Halos 

Ishani Roy^, Chi- Wang Shu^, and Li-Zhi Fang^ 
ABSTRACT 

With a state-of-the-art numerical method for solving the integral-differential 
equation of radiative transfer, we investigate the flux of the Lya photon z/q emer- 
gent from an optically thick halo containing a central light source. Our focus 
is on the time-dependent effects of the resonant scattering. We first show that 



rj I the frequency distribution of photons in the halo are quickly approaching to a 



locally thermalized state around the resonant frequency, even when the mean 



Qh' intensity of the radiation is highly time-dependent. Since initial conditions are 

O ■ forgotten during the thermalization, some features of the flux, such as the two 

^ • peak structure of its profile, actually are independent of the intrinsic width and 

^ ■ time behavior of the central source, if the emergent photons are mainly from 

photons in the thermalized state. In this case, the difference |i^± — t'ol? where u^ 
^ ' are the frequencies of the two peaks of the flux, cannot be less than 2 times of 

Doppler broadening. We then study the radiative transfer in the case where the 



m ■ light emitted from the central source is a flash. We calculate the light curves of 

the flux from the halo. It shows that the flux is still a flash. The time duration 
Q ■ of the flash for the flux, however, is independent of the original time duration 

5 ■ of the light source but depends on the optical depth of the halo. Therefore, the 

spatial transfer of resonant photons is a diffusion process, even though it is not 
a purely Brownian diffusion. This property enables an optically thick halo to 
trap and store thermalized photons around z/q for a long time after the cease of 
the central source emission. The photons trapped in the halo can yield delayed 
emission, of which the profile also shows typical two peak structure as that from 
locally thermalized photons. Possible applications of these results are addressed. 

Subject headings: cosmology: theory - intergalactic medium - radiation transfer 
- scattering 
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1. Introduction 

The transfer of Lja radiation in optically thick medium is fundamentally important for 
the understanding of the physics of halos around Lya photon sources or clouds nearby these 
sources. It includes Lya forest, damped Lya system, Lya blob, Lya emitter and fluorescent 
Lya emission of galaxies and quasars, as well as the optical afterglow of gamma ray bursts. 
The profiles of the emission and absorption of the Lja radiation from these sources are 
powerful tools to constrain the mass density, velocity, temperature and the fraction of neutral 
hydrogen of IGM at various redshifts (Miralda-Escude & Rees 1998; Miralda-Escude 1998; 
Haiman & Cen 2005; Tasitsiomi 2006; Totani et al. 2006; McQuinn et al. 2007). 

It is well known that the resonant scattering of Lja photons with neutral hydrogen 
atoms has a profound effect on the time, space and frequency dependencies of the transfer 
of Lja photons. An analytical solution of the integro-differential equation of the resonant 
radiative transfer revealed that the resonant scattering leads to a local Boltzmann distribu- 
tion of photons in a small frequency range around the Lja frequency uq (Wouthuysen 1952; 
Field 1958, 1959). The temperature of the Boltzmann distribution is equal to the kinetic 
temperature of the neutral hydrogen atoms. The width of the local Boltzmann distribution 
is increasing with time. This is the so-called Wouthuysen- Field effect, which is important in 
21 cm cosmology (Fang 2009). 

Besides Field's solution, no other time-dependent analytical solutions of the integro- 
differential equation are available. All other analytical solutions (Harrington 1973; Neufeld 
1990; Dijkstra et al. 2006) are based on the time-independent Fokker-Planck equation, 
which is the diffusion approximation of the radiation transfer. These time-independent 
solutions are important, but can only be used to describe the "limiting asymptotic behavior" 
of the radiative transfer (Adams 1972). It gives nothing of the time-development of the 
local Boltzmann distribution. The Monte Carlo numerical method has also been used to 
solve the radiative transfer of Lja photons (Lee 1974, Zheng & Miralda-Escude 2002; Ahn 
2002; Cantalupo et al. 2005; Verhamme et al. 2006). Yet, the time-development of the 
Wouthuysen-Field effect is also absent in this approach. 

Time-independent approximation would be reasonable if the length scale / and time 
scales t of the problem considered satisfy the condition l/t <^ c. In this case, one can take 
c — )■ oo, and then, the time derivative term d/cdt of the radiative transfer equation can 
be dropped. The condition l/t ^ c, however, cannot be satisfied if the size of the neutral 
hydrogen halo is large and the time scale of the Lja photon source is small. For instance, 
the time scale tafter of the optical afterglow of GRBs is of the order of a few tens of hours (e.g. 
Tanvir et al. 2009; Salvaterra et al. 2009), or a few days (Vreeswijk et al. 2004), while the 
column number density of neutral hydrogen is as large as IQ^^"^^ cm~^, and number volume 
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density is about 10^ -10^ cm~^ (Vreeswijk et al. 2004). Therefore, the size of the neutral 
hydrogen halos should be much larger than ctaftcr- It is improper to treat the short-time 
problems with the solutions of "limiting asymptotic behavior". 

Recently, a state-of-the-art numerical solver for the kinetic equations has been devel- 
oped. This solver is based on the weighted essentially non-oscillatory (WENO) scheme 
(Jiang & Shu 1996). It has been developed to solve the Boltzmann equations (Carrillo et 
al. 2003, 2006) and radiative transfer (Qiu et al. 2006, 2007, 2008). This numerical solver 
has successfully passed the tests of conservation of photon number. Field's solutions and the 
Wouthuysen-Field effect etc., and has been properly used as the solver of the transfer of 
resonant photons (Roy et al. 2009a, b, c). 

We will study, in this paper, the time- dependent behavior of the Lya radiation transfer 
in an optically thick medium. We will not try to explain any specific observed Lya spectrum, 
instead we will study the physical features of the resonant photon transfer, which can not be 
addressed with previous solvers. For instance, we will show that the frequency distribution 
of the Lya photons can keep in a locally thermalized state even when the intensity and flux 
are highly time- dependent. This feature is essentially important, as thermalization generally 
will lead to the initial conditions being forgotten. Our solver is also able to study the transfer 
of a light flash in optically thick halo. It shows that the nature of the transfer is a diffusion 
process, though it is not a purely Brownian diffusion. This property leads to the trap and 
store of photons thermalized around the Lya frequency for a long time after the cease of the 
central source emission. 

This paper is organized as follows. Section 2 presents the theoretical background of 
the Lya photon transfer in an optically thick medium. Section 3 gives the solution of Lya 
photons emergent from an optically thick spherical halos with a steady source located at 
the center of the halo. Section 4 is on the solutions when the central source is a flash. A 
discussion and conclusion are given in Section 5. The details of the numerical implementation 
is given in the Appendix. 



2. Theory of Ly-alpha radiative transfers in optically thick halos 

2.1. Optically thick halos 

The property of the halo around individual luminous object depends on the luminosity, 
the spectrum of UV photon emission, and the time-evolution of the center object. For 
luminous objects at high redshift, the halos generally consist of three spheres (Cen 2006; 
Liu et al. 2007). The most inner region is the highly ionized Stromgren sphere, or the HII 
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region, which is optically thin of Lja photons. The second region, which is just outside the 
HIT region, is optically thick of Lya photons. The temperature of the baryon gas is about 
10"^ K, which is due to the heating of the UV photons. The third region is outside of the 
heated region. It is un-heated, and therefore, the temperature of the baryon gas can be as 
low as 10^ K. 

In this context, we will study the transfer of Lya photons in a radius R spherical 
halo of neutral hydrogen with temperature T in the range of 10^ to 10^ K. Assuming the 
uniformly distributed HI gas has number density n-m, the optical depth over a light path 
dl is dri^u) = a{i')n}iidl, where (t(z/) is the cross section of the resonant scattering of Lya 
photons by hydrogen, given as 

a{x) = cro(p{x, a) (1) 

where x is the dimensionless frequency defined hj x = {u — z/o)/Az/£), uq being the resonant 
frequency. Auu = vq{vt/c) is the Doppler broadening, and vt = \/2kBT /m. Therefore, x 
measures the frequency deviation Az/ = \v — vq\ in units of Az//) = 1.06 x 10^^(T/10^)^/^ Hz. 
o"o = T^e^ f /iTieC^i'D = 1-10 X 10^^ cm^ is the cross section of the resonant scattering at the 
frequency uq = 2.46 x 10^^ s~^. The function 0(x, a) in equation (1) is the normalized profile 
given by the Voigt function as (Hummer 1965) 

The parameter a in equation (2) is the ratio of the natural to the Doppler broadening. For 
the Lya line, a = 4.7 x 10^^(T/10^)^^'^. The optical depth of the halo with column number 
density of neutral hydrogen Nm = rimR is 

t{x) = Nui<j{x) = To(f){x,a). (3) 

Since 0(0, a) = I/a/tF when a ^ 1, the hne-center optical depth is then r(0) = to/a/tF, and 



.„ = 1.04 X lOM ^V"7-^ I ^ (4) 



2.2. Radiative transfer equation in spherical halo 

The radiative transfer of Lya photons in spherical halo is described by the equation of 
the specific intensity /(?7,r, x,/i) as 

dl dl {l-fi^)dl _ dl _ 
dr] dr r dfi dx 

—(t){x]a)I+ 7l{x,x';a)I{ri,r,x',fi')dx'dfi' + S, (5) 
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where we use dimensionless time 77 defined as rj = cnmcrQt and dimensionless coordinate r 
defined as r = nuiaQVp, with Vp being the physical radial coordinate. That is, rj and r are, 
respectively, in the units of mean free flight-time and mean free path of photon uq. For a 
signal propagated in the radial direction with the speed of light, we have r = rj + const. In 
equation (5) fi = cos 6 is the direction relative to the radial vector r. 

The re-distribution function 7^(x, x'; a) gives the probability of a photon absorbed at the 
frequency x', and re-emitted at the frequency x. It depends on the details of the scattering 
(Henyey 1941; Hummer 1962; Hummer 1969). If we consider coherent scattering without 
recoil, the re-distribution function with the Voigt profile equation (2) is 

lZ{x,x'; a) = (6) 

1 



7,3/2 



e 

\x-x'\/2 



. — 1 / •^min ~r "U \ — 1 / -^inax ^ 

tan — tan 



du 



where Xmin = min(a;, x') and Xmax = max(a;, a;'). In the case of a = 0, i.e. considering only 
the Doppler broadening, the re-distribution function is 

1 
TZ{x,x') = -erfc[max(|x|, |a;'|)]. (7) 

The re-distribution function of equation (7) is normalized as J_ TZ{x,x')dx' = (f){x,0) = 
7i~^^'^e~^ . With this normalization, the total number of photons is conserved in the evolution 
described by equation (5). That is, the destruction processes of Lya photons, such as the 
two-photon process (Spitzer & Greenstein 1951; Osterbrock 1962), is ignored in equation 
(5). In equations (6) and (7), we also do not consider the recoil of atoms. It is equal to 
assume the mass of atom is very large. The effect of recoil actually is under control (Roy et 
al. 2009c). We will address it in next section. 

In equation (5), the term with the parameter 7 is due to the expansion of the universe. 
If Tin is equal to the mean of the number density of cosmic hydrogen, we have 7 = r^^p, and 
Top is the Gunn- Peterson optical depth. Since Gunn- Peterson optical depth is of the order 
of 10® at high redshift (e.g. Roy et al. 2009c), we will take the parameter 7 = 10^^ — 10^®. 



2.3. Eddington approximation 

When the optical depth is large, we can take the Eddington approximation as 

I{ri, r, X, fj,) ~ J(?7, r, x) + ?>jiF{rj, r, x) 



1 r+i 



where J{rj, r, x) = | /_^ /(r^, r, x, iJi)dji is the angularly averaged specific intensity and -F(?7, r. 



X] 



1 r+i 



2 



/-^ fil{ri,r,x, fi)dfj, is the flux. Defining j = r^J and / = r^F, Eq.(5) yields the equations 
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of j and / as 

Ol O f f Ol 

— + — = -0(a;; a)j + / 7^(a;, x'; a)jdx' + 7— + r^S, (9) 

The mean intensity j{r],r,x) describes the x photons trapped in the halo by the resonant 
scattering, while the flux f{ri,r,x) describes the photons in transit. 

The source term S in the equations (9) and (10) can be described by a boundary 
condition of j and / at r = tq. We can take tq = 0, as tq is not important if the optical 
depth of the halo is large. Thus, we have 

j{r], 0, x) = 0, f{r], 0, x) = S'o0s(x), (11) 

where Sq, and 0s(x) are, respectively, the intensity and normalized frequency profile of the 
sources. Since equation (9)- (11) are linear, the intensity 5*0 can be taken as any constant. 
That is, the solution f{x) of Sq = S is equal to Sfi{x), where fi{x) is the solution of S* = 1. 
On the other hand, the equations (9) and (10) are not linear with respect to the function 
(ps{x), i.e. the solution f{x) of 0s(x) is not equal to 0s(x)/i(x), where /i is the solution of 
(ps{x) = 1. 

In the range outside the halo, r > R, no photons propagate in the direction fi < 0. 
Therefore, the boundary condition at r = R given by /^ fiJif], R, x, fi)dfi = is then (Unno 
1955) 

j{r],R,x)=2f{r],R,x). (12) 

There is no photon in the field before t = 0. Therefore, the initial condition is 

j(0,r,x) = /(0,r,x) = 0. (13) 

We solve equations (9) and (10) with the numerical method developed recently (Roy et 
al. 2009a, 2009b, 2009c). Some details of this method is given in the Appendix. We first 
solve the problems when the sources S is steady. 



3. Solutions of steady sources 

3.1. Time scale of escape 

First we consider steady sources. That is, the parameter So in eq.(ll) is time-independent 
after the switch-on of the sources at r^ = [eq.(13)]. A typical solution of the flux f{ri,r,x) 
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given by equations (9) and (10) is shown in Figure 1, for which the source is taken to be 
5*0 = 1 and 0s(a;) = (l/y^)e~^ , i.e. the emission hne width is equal to the Doppler broad- 
ening. The parameters a and 7 are taken to be 10~^ and 10~^, respectively. The effect of 
7 = 10~^ actually is ignorable in these solutions. The left panel is the solutions /(?7,r, x) 
at radius r = 10^ and time 77 = 500, 1000, 2000 and 3000 with the boundary condition 
equations (11) and (12). The solutions approach to a stable state at time 77 > 2000. 





Fig. 1. — Left panel is the solution of flux /(r^, r, x) of equations (9) and (10) at r = i? = 10^ 
and 7] = 500, 1000, 2000 and 3000 with the boundary condition equation (12), R = 10^. The 
source is 6*0 = 1 and (f)s{x) = {l/^/^T)e~^ . The parameters a and 7 are taken to be a = 10^^ 
and 7 = 0^^. The right panel is the time-dependence of the total flux Ft{ri,r) at r = 10^. 

The right panel of Figure 1 is the total flux Ft{r]) = J f{ri, r, x)dx. It shown again that 
total flux approaches to a stable state with F^ = 1 at time t] > 2000. If the Lya photon 
transfer is due only to spatial diffusion, the time scale of a spatial transfer with an optical 
depth To = 10^ should be as large as 77 ~ Tq = 10"^. However, Figure 1 shows that photons 
have already escaped from the tq = 10^ hole within the time rj ~ 2000, which is much less 
than the time scale of a purely Brownian diffusion. Therefore, the transfer of photons should 
not be a process of purely Brownian diffusion. This point is very well known in early studies 
on the escape of resonant photon from opaque clouds (Osterbrock 1962; Harrington 1973; 
Avery & House 1968). The time scale of escape shown in Figure 1 is also consistent with the 
estimation given by Monte Carlo simulations (Adams, 1975; Bonilha et al. 1979). However, 
these works are mainly based on the time scale of the escape of photons with frequency, at 
which the photon can take a "single longest excursion" (Adams 1972). Figure 1 shows that 
the escape time scale 77 = 2000 is available not only for photons which can take a "single 
longest excursion", but also for all photons with frequency around z/q or a; = 0. 



The right panel of Figure 1 also plays the role of testing our algorithm. Because 
/ (f){x, a)jdx = f R{x, x'\ a)jdxdx' , eq.(9) yields 

Fi(77,r)=Fi(r7,0) = 5o (14) 

when the solution approached the stable state. Equation (14) is the conservation of photon 
number. Since numerical errors which have accumulated over a long time evolution could be 
huge, eq.(14) is useful to check the reliability of the code. The right panel of Figure 1 shows 
perfectly Ft{ri,r = 10^) = Ft{ri,0) = 5*0 = 1 at stable state. 



3.2. Time scale of local thermalization 



Figure 2 shows a solution of the mean intensity j(?7,r, x) of equations (9) and (10) at 
r = 10^ and time rj = 200, 300, 500. The source is the same as Figure 1, namely 5*0 = 1 and 
0s(x) = {l/^/^T)e~^ . Other parameters are also the same as Figure 1. 
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Fig. 2. — Mean intensity j{r], r, x) at r 
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be a = 10 ^ and 7 = 0^ 



10^ of equations (9) and (10) at time r] = 200, 300, 
^\ -o- rjj_^g parameters a and 7 are taken to 



A remarkable feature of the solutions is to show a flat plateau around x = 0. As has 
been shown by Roy et al. (2009c), the fiat plateau actually is the Boltzmann statistical 
equilibrium distribution around x = when the atomic mass is infinite. If the mass is finite, 
i.e. considering the recoil in the re-distribution functions (6) or (7), the fiat plateau will 
become e~^^^, where b = huQ/mvTC. This is the local Boltzmann distribution required by 
the Wouthuysen-Field effect (Roy et al. 2009b). 
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Without resonant scattering, the pure absorption will lead to the flux of z/q photons at 
r = 10^ to be attenuated by a factor e^*-^-' ^ 10~^^. Therefore, Figure 2 shows that a major 
effect of resonant scattering is to restore uq photons in optically thick halos. According to 
the re-distribution function TZ{x, x') equation (5), the probability of transferring a x' photon 
to a I a; I < \x'\ photon is larger than that from x' to |x| > \x'\. Therefore, the net effect of 
resonant scattering is to bring photons with frequency x' 7^ to the central range a; ~ 
of frequency space. Photons of x ~ are then effectively restored. Moreover, the restored 
photons are thermalized. We see from Figure 2 that in the time range from 77 = 200 to 500, 
the mean intensity j at x = quickly increases by a factor larger than 10. In the same 
time, the flat plateau or the local thermalization around x = is perfectly held. Therefore, 
the time scales tther of uq photon restoring and thermalization should be less than t] = 200, 
which is much smaller than the time scale t escape of escaping from a halo of r = 10^. 
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Fig. 3. — Flux f{ri,r,x) at r = 10^ and time rj = 200, 300 and 500, with sources left panel: 
5*0 = 1 and 0s (x) = (l/v/vr)e~^ , and right panel: 5*0 = 1 and 0s(x) = {l/>y7T/2)e~'^^ . 
Other parameters are the same as in Figure 1 



The result of tther ^ t escape i^ optically thick halos is very important. One can conclude 
that the photons of the flux / around x = emergent from an optically thick halos actually 
come from the thermalized photons, regardless of the initial distribution of the photons. The 
initial conditions of the photon source have already been forgotten during the thermalization. 
Therefore, one can expect that the proflle of flux / has to be independent of the sources. 
We demonstrate this point with Figure 3. 



Figure 3 presents the flux / given by equations (9) and (10) with the same parameters 
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as the solution of Figure 1, and the source profiles are (j)s{x) = (l/-\/7r)e~^ (left) and (f)s{x) = 
(l/A/7r/2)e~^^ (right). That is, the line widths are, respectively, 1 and \j\pl. We see that 
the left and right profiles of Figure 3 are almost identical within |x| < 3. Therefore, the 
initial distribution of photon frequency is already forgotten, and the photons of the left and 
right cases at r^ > 200 actually are from the same thermalized sources. The profiles of the 
flux are also held if the line width is broader than 1. We will show this point when the source 
has a continuant spectrum (§3.4). 



3.3. Two peaks in the flux profile 

The profiles of either j or / shown in Figures 1-3 have two peak structure. The flux 
/ is dominated by photons with frequency x± ~ ±(2 — 3). The two peak structure has 
been recognized in all the time-independent solutions of the Fokker-Planck approximation 
(Harrington 1973; Neufeld 1990; Dijkstra et al. 2006), and Monte Carlo simulations (Lee 
1974; Zheng & Miralda-Escude 2002; Ahn et al. 2002; Cantalupo et al. 2005; Verhamme et 
al. 2006). A point we would like to emphasize is that this structure is independent of the 
profile of the source S. It is because the initial properties of the central sources have been 
forgotten during the local thermalization. 

Since the shape of the locally thermalized distribution is time- independent, the fre- 
quency of the two peaks, |a;±|, at a given r is also time-independent. When the photons 
of the flux / mainly come from the locally thermalized photons, the frequency of the two 
peaks, |x±|, should not be described by the relation |a;±| = {(itY^^ , because this relation is 
from a solution of the time-independent Fokker-Planck equation, which does not describe 
the thermalization (Neufeld 1990). 

Figure 4 presents the peak frequency |x±| as a function of ar for solutions given by 
equations (9) and (10) with a = 10~^. We consider only r > 10^, as in the case r < 
10^ photons do not undergo a large enough number of scattering, and therefore, are not 
completely thermalized yet. With the dimensionless variables, ar is equal to ar. A line 
\x±\ = {arY^^ is also shown in Figure 4. It shows that our numerical result of \x±\ is 
significantly different from the (aro)^/^-law if ar < 30. That is, in the range ar < 30, 
photons of the fiux / are dominated by the locally thermalized photons. The frequency x± 
actually is dependent on the width of the flat plateau or the locally thermalized distribution 
of j. Therefore, x± is always larger than two. This feature has also been addressed in Bonilha 
et al (1979). The curve of Figure 4 is approximately available for a = 10~^ and 10~^. Thus, 
|a;±| is larger than {ary^^ if r < 3 x 10^ and 3 x 10^ for a = 10~^ and 10~^, respectively. 



- 11 - 




Fig. 4. — The positions of the peaks |a;+| as a function of (arY'^ for the solutions of eqs. (9) 
and (10) with a = 10~^. Other parameters are the same as in Figure 1. The dashed hne is 

for x± = ±(ar)^/^. 



Figure 4 shows a slowly increase of |a;-|-| with r in the range ar < 30, and then, it 
approaches {arY^^ for larger ar. When r is large, more photons of the flux are attributed to 
the resonant scattering by the Lorentzian wing of the Voigt function. |x-|- 1 is then approaching 
to (arY^^- It should be emphasized once again that all these results are independent of the 
intrinsic width of Lya emission from the central source. 



3.4. Absorption spectrum 

If the radiation from the sources has a continuant spectrum, the effect of neutral hydro- 
gen halos is to produce an absorption line at z/ = z/q- The profile of the absorption line can 
also be found by solving equations (9) and (10), but replacing the boundary equation (11) 
by 

j(r/,0,x) = 0, f{i],0,x) = So. (15) 

That is, we assume that the original spectrum is fiat in the frequency space. 

A solution of the time evolution of j and / at r = 10^ with the source equation (15) 
is shown in Figure 5. The optical depths at the frequency \x\ > 4 are small, and therefore, 
the Eddington approximation would no longer be proper. However, those photons do not 
strongly involve the resonant scattering, and hence they do not significantly affect the so- 
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Fig. 5. — Solutions of j{ri,r,x) (left) and f{ri,r,x) (right) of eqs. (9) and (10) at r = 10^ 
and 7] = 200, 300 and 500. The source is given by eq. (16). The parameters are a = 10^^ 
and 7 = 10~^. 



lution around x = 0. Therefore, the solution is still useful to study the profiles of j and 
/ around a; = 0. The profile of / typically is an absorption line with width given by the 
position of the two peak structure. As expected, the profile in the range \x\ < 4 is the same 
as the left panel of Figure 1. It shows again that the two peak structure is independent of 
the frequency spectrum of the central source. 

The flux / of Figure 5 has a fiat wing in both sides of x > 3 and x < —3, because the 
effect of resonant scattering is negligible for photons with frequency |x| > 4, and they can 
freely transfer from r = to 10^. On the other hand, the mean intensity j, which is the 
isotropic component of the intensity J [eq.(8)], does not have wings at |x| > 4. That is, 
resonant scattering cannot store photons with |a;| > 4 in the halo. Nevertheless, the mean 
intensity j still has a flat plateau around x = 0. It means that in the period from rj = 200 
to 1] = 500, the halo trapped and stored more and more photons of |x| < 4. These photons 
are in the locally thermalized state. 



3.5. The estimation of the HI column density 

As an application of the absorption spectrum in Figure 5, we study the profile of the red 
damping wing of the optically thick IGM at high redshifts. Since the variable x is independent 
of redshift, the profile of a red damping wing is directly given by the flux f{ri,r,x) at the 
wing of low frequency a; < 0. 
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Fig. 6. — Red damping wing of (a) the DLA model eq.(16) tq = 2 x 10^ (dashed hne); (b) 
the solution fix) of §3.4 at r = 10^ and r^ = 9 x 10^ (solid line). For both (a) and (b), the 
parameter a = 10^^. 



If the hydrogen clouds are located far from the Lya sources, the red damping wing can 
be modeled as the absorption of an optically thick halos. The red damping wing of damped 



Lja system (DLA) model is then given by f{x) 
Voigt function equation (2) as 



e '^^^' and x < 0, where t{x) is from the 



t[x) = rQ(f){x, a) 



To- 



TT 



3/2 



dy 



(y 



x)'^ + a? 



(16) 



The column number density of neutral hydrogen atoms, A^m? generally is estimated with 
fitting profile equation (16) with observation, and then, A^hi can be found from tq by equation 
(4). If the light source is located in a hydrogen cloud, the column number density given by 
the fitting of equation (16) should be underestimated, because resonant scattering helps the 
transfer of resonant photons. 

As an example Figure 6 gives (a) the red damping wings of eq.(16) with tq = 2 x 10^, 
and (6) a solution / of §3.4 at r = 10^, or tq = 10^ and rj = 9 x 10^. The profiles of (a) and 
(6) are very different from each other. The former is quickly dropping when |a:| is less than 
about 3, while the later at |x| < 3 is softly dependent on x. The curve of (6) approaches to 
a much higher bottom at a; = 0. 

More importantly. Figure 6 shows that the curve (a) with tq = 2 x 10^ is more or less 
close to the curve (6). That is, the DLA model at tq = 2 x 10^ may give a similar data fitting 
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as the solution with resonant scattering. Therefore, for optical depth tq = 10^, the DLA 
model of equation (16) will underestimate the column number density of neutral hydrogen 
by about two orders. 



4. Solutions of flash sources 

4.1. Frequency-dependence of Ly-alpha transfer 

If the light source is significantly time-dependent, like the optical afterglow of GRBs, 
we can treat the source as a flash described by a boundary conditions as: 

,(,,0,.) = 0, /(„,0,.) = /^»^-'^'- »<"<"» (17) 

10, r]> 7]Q. 

It describes a flash within a time range < r^ < t^Qj or the time duration is A?] = tjq. We 
consider the case of r > tjq. That is, the size of the halo is larger than the spatial range 
lasted by the flash, as photon spatial transfer in optically thick medium cannot be faster 
than the speed of light. The initial condition is still the same equation (13). 

Figure 7 presents two time-dependent solutions of the mean intensity j and flux /: one 
is for a flash of equation (17) with r^o = 100 at r = 10^; the other is for a flash with r^o = 500 
at r = 10'^. Other parameters are 5*0 = 1, a = 10^^ and 0s(a;) = (l/y^)e^^ . We still see 
the typical flat plateau of j in all times, even when the original time duration of the flash is 
as short as Ar] = 100. 

The time dependence of j has a rising phase and a decaying phase. The thermalization 
of j is held in both rising and decaying phases. We see from Figure 7 that the rising and 
decaying phases are frequency-dependent. Photons at the two peaks reach their maximum 
at about rj = 200 (top right panel) and r] = 4000 (bottom right panel), while the flat plateau 
reaches their maximum at about r] = 500 (top right panel) and t] = 6000 (bottom right 
panel). That is, the halo holds a locally thermalized photons for a much longer time than 
the original time durations At] = 100 (top) and Arj = 500 (bottom). 

The time-evolution of / also consists of a rising phase and a decaying phase. Therefore, 
the flux emergent from the halo is also a flash. However, the time duration is very different 
from the original one. For the top panel, we see that the proflle of / is almost time- 
independent from f] = 300 to 500. That is, the time duration 500 — 300 = 200 is much 
larger than the original one At] = 100. For the bottom panel, the time-independent flux is 
from 1] = 3000 to 5000, and therefore, the time duration of the flash is about 2000, which is 
also much larger than the original time duration rj = 500. 
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Fig. 7. — The profiles oi j{ri,r,x) and f{ri,r,x) with time- dependent source [eq.(17)]. Top 
paneh j{ri,r,x) (left) and f{ri,r,x) (right) of r]Q = 100 at r = 10^. Bottom panel: r]Q = 500 
at r = 10^. The parameters a = 10~^ and 7 = 10~^. 



4.2. The light curves of a flash 



To further study the feature of time dependence of flash, we plot Figure 8, which gives 
the light curve of the flux / at a; = (top panels), and the total flux (bottom panels). The 
top panels show the light curve of rising and decaying phases of the z/q photon flux /. With 
these light curves, one can find the maximum of the flux / at time rjmax] the rising phase is 
then T] < rjmax, and decaying phase is 77 > rjmax- For each curve, one can also define a time 
duration A77 = ?72 — "'^i, where % < rjmax and ri2 > rjmax, and both are given by the condition 
/(^i,2, r, 0) = 0.9 f{rimax, r, 0). The time duration is then Ar] = rj2 — tji. 

The two top panels of Figure 8 are for a flash source with original time duration A?] = 50. 
Figure 8 shows that the time duration of the flash. A?], is r-dependent. At r = 0, i.e. at the 
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Fig. 8. — Top panel: the light curves of f{ri,r,x) at x = at r = 10^ (left) and r = 10'^ 
(right) for flash source with rjQ = 50. Bottom panel: the light curves of the total flux 
Ft{r],r) = J f{ri,r,x)dx at r = 10^ (left) and r = 10^ (right). 



source, Arj is 50. At r = 10^ (top left of Figure 8), Arj is about 200, while at r = 10^, we 
have At] c^ 2000. Therefore, the time duration Arj is increasing with r. This result shows 
that the time duration seems to be dependent mainly on the size r (or optical depth r) of 
the halo, regardless the initial time duration Arj = rjo- 

The bottom two panels of Figure 8 are the light curves of the total flux Ft{ri,r) = 
J f{r],r,x)dx. The peak times r]max of total flux generally are less than that of uq photon 
flux, as the z/q photon needs longer restoration time. The time durations given by the total 
flux are also a little less than that of z/q photon flux, but it still shows that the time duration 
is mainly dependent on the size r of the halo. 

Figure 9 presents the light curves of flash sources with time duration A77 = 770 = 1, 
5 and 20. The halo size is r = 10^. For the case of rjQ = 1, we have rjQ <^ r. Therefore, 
the source can be considered as a pulse. The top panels of Figure 9 are for flux of z/q 
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Fig. 9. — Top panel: the time-dependence of f{r],r,x) at a; = and at r = 10^ of flash 
source with rjo = 1 (left), 5 (middle) and 20 (right). Bottom panel: the light curves of the 
total flux Ft{ri,r) of the corresponded top panel. 



photons, while the bottom panels are the corresponded total flux. Although the three flash 
sources have very different time durations at r = 0, their light curves of / at x = are very 
similar. The maximum values of / for r^o = 1, 5 and 20 can even be described by relations 
as /20 — 4/5 ^ 20/1. The coefficients 4 and 20 are from the ratio of the total numbers of 
photons of the three flashes. The three light curves of Ft at the bottom panels of Figure 9 
are also similar from each other, although they are not as good as the top three curves /. 
This is because the light curves are frequency-dependent. Nevertheless, we still see the three 
maximum values of Ft also roughly satisfy the relation F20 — 4F5 ~ 20-Fi. 

Either Figure 8 or Figure 9 reveals that the time scales of the propagation of a flash in 
halos are mainly dependent on the size r, regardless the original time duration. This feature 
indicates that the spatial transfer of a flash essentially is a diffusion process. As mentioned 
in §3, the spatial transfer of resonant photons cannot be described as a purely Brownian 
diffusion process, by which the time duration Arj should increase with r^ or r^. On the 
other hand, if the spatial transfer of a photon can be, in average, described by constant 
speed, the time duration Arj of a flash should be, in average, equal to the original time 
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duration, or at least, dependent on the original time duration. However, Figures 8 and 9 
shows that the time duration At] is roughly proportional to r, and independent of the initial 
time duration. Therefore, the spatial transfer of resonant photons essentially should still be 
a diffusion process. 



4.3. Delayed emission line of stored photons 



The effect of trapping and storing photons in an optically thick halo can be more clearly 
revealed with a flash source. Let us consider a flash source with a continuant spectrum, which 
is described by equation (17), but we take S^cpsix) = 1. The time evolutions of j and / for 
r^o = 50 at r = 10^ are presented in Figure 10, in which we take the time to be r^ = 200, 300 
and 500. 









Fig. 10. — Profiles oi j{r],r,x) (top) and f{ri,r,x) (bottom) at r = 10^ when a fiash source 
eq.(17) with rjQ = 50 and SQ(f)s{x) = 1. The time is rj = 200 (left), 300 (middle) and 500 
(right). 



The top panels show the evolution of the mean intensity j in the halo. In an early time 
1] = 200, there are photons in the central range |x| < 4 as well as in wings |a;| > 4. At later 
time f] = 300, wing photons disappear, because all wing photons from the fiash source have 
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already escaped from the r halo. At time rj = 500, the mean intensity j of |a;| < 4 is still 
about the same as j at rj = 300. The flat plateau of j is shown in all time. Therefore, the 
locally thermalized photons are stored in the halo at least from time rj = 200 to 500, which 
is much longer than the r^o = 50. 

The evolution of the flux / given by the bottom panels of Figure 10 is more interesting. 
At the time r] = 200, the flux shows a typical absorption line at z/q. However, at r] = 300, 
the flux / becomes a typical emission line with two peak profile. At time r] = 500, / is still a 
two peak emission line. The flux of the emission ai rj = 500 is as strong as that at 77 = 300. 
Its light curve is similar to that of Figures 8 and 9. Note that Figures 8 and 9 are for a 
source of emission line, while Figure 10 is for a continuant spectrum. The similarity of the 
light curves of Figure 10 with Figures 8 and 9 is again due to the local thermalization and 
the diffusion in the physical space, both of which lead to the initial frequency spectrum and 
the time dependence of the photon sources being forgotten. 

The emission aX rj > 300 is a delayed emission, as the flash of source has already ceased. 
The photons of the delayed emission is provided by the |x| < 4 photons stored in the halo 
r < 10^. The time duration of the delayed emission is about the same as the time scale of 
the decaying phase of Figures 8 and 9. Therefore, it is also proportional to r, regardless the 
original time duration. Thus, at large r, a flash with a continuant spectrum and very short 
original time duration rjo can produce a two peak emission with time duration proportional 
to r. 



5. Discussions and conclusions 

The resonant scattering made the transfers of resonant photons in physical space and 
frequency space to be coupled from each other. It leads to the time scale r] of the spatial 
transfer of the resonant photons in the halo with optical depth r ^ 1 being much faster 
than a purely Brownian diffusion process requiring 77 oc r^. However, essentially it is still a 
diffusion process, which can be approximately described by r^ oc r^^^, with the index H less 
than but very close to 1. It is possible, if we consider the single longest excursion playing 
the role of a long-range dependence, that this diffusion process has a positive correlation, or 
is a fractal Brownian diffusion (e.g. Beran 1994). 

The number of photons basically is conserved. Thus, an optically thick halo is a store 
of photons with frequency ~ z/q. The time scale of the store is the same as the time scale 
of the above-mentioned diffusion, i.e. approximately proportional to the optical depth r. 
Moreover, the stored photons are always in the state of local Boltzmann distribution, even 
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when the mean intensity is highly time-dependent. The initial conditions are forgotten in the 
process of approaching the locally thermal equilibrium. The local Boltzmann distribution is 
independent of the frequency spectrum and time-dependence of the source. 

All these features show that the z/ ~ z/q photons play the central role of radiation transfer 
with resonant scattering. The major difference between our solutions and some analytical 
solutions (Harrington 1973; Neufeld 1990; Dijkstra et al. 2006) is also around a; = 0. Since 
the analytical solutions are based on the assumption (j){x) = a/y^x^, and the Gaussian core 
g-x jg ignored, it generally leads to J{x = 0) = 0. With these approximations the solutions 
cannot show the effects of restoration and thermalization of photons around x = 0. 

These basic properties found in our numerical solutions yield the following features of 
Lya photons emergent from an optically thick halo. 

1. At a given r, the profile of the two peaks of the fiux is time-independent, and also 
independent of the initial profile of the photons. Therefore, it is impossible to estimate 
the line profile of the source from the profile of the flux emergent from an optically 
thick halo. 

2. The frequencies |a;±| of the two peaks of the flux is not less than about 2. This would 
be useful to estimate the kinetic temperature of the neutral hydrogen atoms. 

3. The resonant scattering makes the flux of the red damping wing very different from 
that of the DLA model. The flux is non-zero at frequency z/q or a; = 0. These results 
would be useful to discriminate the DLA model with models which consider the effect 
of resonant scattering. 

4. The time scales of light curves of the delayed emission of Lya photons from a flash 
source is mainly determined by the optical depth of the halos. On the other hand, the 
halo is transparent for high energy photons. A comparison between the light curves of 
Lya photons and high energy photons would be useful to detect the halo. 

For halos with large optical depth, the parameter 7 is small even at high redshift. 
When 7 ~ 10~^, the effect of cosmic expansion on photon evolution in the frequency space 
actually is negligible. All solutions f{ri,r,x) given in this paper are almost independent of 
the Hubble expansion. The effect of 7 would be important if the considered time scale rj of 
/ is comparable with that of the Hubble expansion. 

This work is supported in part by NSF grant AST-0506734 and ARO grant W911NF- 
08-1-0520. 
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A. Numerical algorithm 

To solve equations (9) and (10) as a system, our computational domain is {r,x) e 
[0, rmax] X [xieft, bright], where Vmax, Xieft and Xright are chosen such that the solution vanishes 
to zero outside the boundaries. We choose mesh sizes with grid refinement tests to ensure 
proper numerical resolution. In the following, we describe numerical techniques involved in 
our algorithm, including approximations to the spatial derivatives, integrals in the frequency 
domain, numerical boundary condition and time evolution. 



A.l. The WENO algorithm: approximations to the spatial derivatives 

The spatial derivative terms in equations (9) and (10) are approximated by a fifth order 
finite difference WENO scheme. 

We first give the WENO reconstruction procedure in approximating g|, 

dj{ri'',ri,Xj) 1 



dx Ax 



{hj+i/2-hj^i/2), (Al) 



with fixed r] = r]"' and r = ri. The numerical fiux /ij+1/2 is obtained by the fifth order WENO 
approximation in an upwind fashion, because the wind direction is fixed. Denote 

h, = jir, n, X,), J = -2, -1, ■ ■ ■ , AT + 3 (A2) 

with fixed n and i. The numerical fiux from the WENO procedure is obtained by 

hj+1/2 = ^lhfli/2 + ^2/if+i/2 + W3/if+i/2, (A3) 

where h^{/2 are the three third order fiuxes on three different stencils given by 

hj+i/2 = ~g^i-i + g^i + 3^i+i' 
"i+1/2 - 3'^i + 6 ■''^^ ~ g'^i+2' 



;(3) 



11. 7, 1 



hj+i/2 - Y^i+i ~ Q^^+'^ "^ 3^i+3) 



and the nonlinear weights Um are given by. 



Um = , Ul = J-—^, (A4) 



^ (^ + A 

1=1 
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where e is a parameter to avoid the denominator to become zero and is taken as e = 10 ^. 
The hnear weights 7; are given by 

3 3 1 

71 = Y^, 72=5, 73 = ^, (A5) 

and the smoothness indicators /?/ are given by, 

13 1 

13 1 

/32 = —{h,~2h,+^ + hj+2f + -{h,-h,+2)\ 

13 1 

To approximate the r-derivatives in the system of equations (9) and (10), we need to 
perform the WENO procedure based on a characteristic decomposition. We write the left 
hand side of equations (9) and (10) as 

Mt + Avir (A6) 

where u = {j, f)^ and 



A 



1 

1 



is a constant matrix. To perform the characteristic decomposition, we first compute the 
eigenvalues, the right eigenvectors, and the left eigenvectors of A and denote them by. A, 
R and R~^. We then project u to the local characteristic fields v with v = R'^u. Now 
ut + Aur of the original system is decoupled as two independent equations as vt + Av^. 
We approximate the derivative v^ component by component, each with the correct upwind 
direction, with the WENO reconstruction procedure similar to the procedure described above 
for o^. In the end, we transform v^ back to the physical space by u^ = Rvr- We refer the 
readers to Cockburn et al. 1998 for more implementation details. 



A. 2. Adaptive mesh procedure for non-uniform grid 

A fifth order conservative finite difference WENO scheme can only be applied to a 
uniform grid or a smoothly varying grid. A smooth transformation. 



e = m (A7) 
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gives us a uniform grid in a new variable C,- In this case ^ is sufficiently smooth, i.e., it has 
as many derivatives as the order of accuracy of the scheme. Therefore the left hand side of 
the (9) and (10) as 

Ut + Aur (A8) 

where u = (j, f)'^ and 

■ 1 



is transformed to 

Ut + A^rU^ (A9) 

and the WENO derivative approximation is now applied to u^. 
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